path_proj = here::here()
path_source = file.path(path_proj, "source")
# source
source(file.path(path_source, "simulation", "simulations_functions_final.R"))
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:reshape2':
##
## smiths
## The following object is masked from 'package:rstan':
##
## extract
# source from other paper
source(file.path(path_source, "funcs_german_paper", "plotReportingTriangle.R"))
# data import
SARI <- as.data.frame(read.csv(file.path(path_proj, "data", "raw", "clean_data_srag_epiweek_delay_table_PR.csv"),))
# set week
set_week_start("Sunday")
SARI$date <- get_date(week = SARI$epiweek, year = SARI$epiyear)
fixed_q <- file.path(path_proj, "source", "models",
"trunc", "1.stan_model_fixed_q_trunc.stan")
fixed_b <- file.path(path_proj, "source", "models",
"trunc", "2.stan_model_fixed_b_trunc.stan")
b_poly <- file.path(path_proj, "source", "models", "trunc",
"3.stan_model_b_poly_trunc.stan")
b_spline <- file.path(path_proj, "source", "models", "trunc",
"4.stan_model_b_spline_trunc.stan")
compiled_models <- list(
fixed_q = cmdstan_model(fixed_q),
fixed_b = cmdstan_model(fixed_b),
b_poly = cmdstan_model(b_poly),
b_spline = cmdstan_model(b_spline)
)
source(file.path(path_source, "functions", "plot_function.R"))
source(file.path(path_source, "functions", "fit_function.R"))
# setting: delay, number of days
seed <- 123
D <- 15
N_obs <- 200
# try-on in region 41002
data_41002 <- SARI %>% filter(regionalsaude == 41002) %>%
select(-c("epiweek","epiyear","regionalsaude")) %>%
relocate(Notifications, .after = last_col())
rownames(data_41002) <- data_41002$date
#transfer to cumu matrix
data_41002[,c(1:27)] <- cumulative_matrix(as.matrix(data_41002[,c(1:27)]))
# now <- as.Date("2024-02-01")
scoreRange <- seq(as.Date("2009-07-20"),as.Date("2009-08-30"),by="7 day")
# input the true case vector
case_true <- as.matrix(data_41002[, c("Notifications")])
rownames(case_true) <- as.character(data_41002$date)
out <- nowcasting_moving_window(data_41002[, c(1:21)], scoreRange = scoreRange,
case_true = case_true,
start_date = as.Date("2009-07-01"),
#start_date = scoreRange[1] - weeks(20),
D = 20, sigma_b = 0.5, seeds = 123,
models_to_run = c("fixed_q", "fixed_b", "b_poly","b_spline"),
compiled_models = compiled_models,
iter_sampling = 5000, iter_warmup = 1000, refresh = 0,
num_chains = 3, suppress_output = T)
## ====================
## now=2009-07-20 (1/6)
## ====================
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 0.8 seconds.
## Chain 2 finished in 0.9 seconds.
## Chain 3 finished in 0.9 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 0.9 seconds.
## Total execution time: 2.9 seconds.
##
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 0.3 seconds.
## Chain 2 finished in 0.3 seconds.
## Chain 3 finished in 0.3 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 0.3 seconds.
## Total execution time: 1.2 seconds.
##
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 0.6 seconds.
## Chain 2 finished in 0.3 seconds.
## Chain 3 finished in 0.4 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 0.4 seconds.
## Total execution time: 1.6 seconds.
##
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 0.4 seconds.
## Chain 2 finished in 0.5 seconds.
## Chain 3 finished in 0.5 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 0.5 seconds.
## Total execution time: 1.7 seconds.
##
## ====================
## now=2009-07-27 (2/6)
## ====================
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 1.4 seconds.
## Chain 2 finished in 1.2 seconds.
## Chain 3 finished in 1.1 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 1.2 seconds.
## Total execution time: 4.0 seconds.
##
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 0.4 seconds.
## Chain 2 finished in 0.4 seconds.
## Chain 3 finished in 0.6 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 0.4 seconds.
## Total execution time: 1.6 seconds.
##
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 0.5 seconds.
## Chain 2 finished in 0.5 seconds.
## Chain 3 finished in 2.3 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 1.1 seconds.
## Total execution time: 3.5 seconds.
##
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 0.6 seconds.
## Chain 2 finished in 0.7 seconds.
## Chain 3 finished in 0.6 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 0.6 seconds.
## Total execution time: 2.1 seconds.
##
## ====================
## now=2009-08-03 (3/6)
## ====================
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 2.9 seconds.
## Chain 2 finished in 2.5 seconds.
## Chain 3 finished in 2.2 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 2.5 seconds.
## Total execution time: 7.8 seconds.
##
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 1.7 seconds.
## Chain 2 finished in 1.6 seconds.
## Chain 3 finished in 1.4 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 1.6 seconds.
## Total execution time: 5.0 seconds.
##
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 1.9 seconds.
## Chain 2 finished in 2.2 seconds.
## Chain 3 finished in 1.9 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 2.0 seconds.
## Total execution time: 6.4 seconds.
##
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 2.2 seconds.
## Chain 2 finished in 3.0 seconds.
## Chain 3 finished in 2.5 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 2.6 seconds.
## Total execution time: 8.1 seconds.
##
## ====================
## now=2009-08-10 (4/6)
## ====================
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 3.3 seconds.
## Chain 2 finished in 3.3 seconds.
## Chain 3 finished in 3.1 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 3.2 seconds.
## Total execution time: 9.9 seconds.
##
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 3.2 seconds.
## Chain 2 finished in 2.6 seconds.
## Chain 3 finished in 2.6 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 2.8 seconds.
## Total execution time: 8.9 seconds.
##
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 5.2 seconds.
## Chain 2 finished in 4.6 seconds.
## Chain 3 finished in 5.0 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 5.0 seconds.
## Total execution time: 15.3 seconds.
##
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 4.6 seconds.
## Chain 2 finished in 4.4 seconds.
## Chain 3 finished in 4.4 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 4.5 seconds.
## Total execution time: 13.8 seconds.
##
## ====================
## now=2009-08-17 (5/6)
## ====================
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 4.1 seconds.
## Chain 2 finished in 4.2 seconds.
## Chain 3 finished in 4.5 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 4.3 seconds.
## Total execution time: 13.1 seconds.
##
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 3.7 seconds.
## Chain 2 finished in 4.4 seconds.
## Chain 3 finished in 4.2 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 4.1 seconds.
## Total execution time: 12.6 seconds.
##
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 9.5 seconds.
## Chain 2 finished in 9.4 seconds.
## Chain 3 finished in 10.3 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 9.7 seconds.
## Total execution time: 29.6 seconds.
##
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 5.3 seconds.
## Chain 2 finished in 6.0 seconds.
## Chain 3 finished in 7.8 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 6.4 seconds.
## Total execution time: 19.3 seconds.
##
## ====================
## now=2009-08-24 (6/6)
## ====================
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 4.5 seconds.
## Chain 2 finished in 4.5 seconds.
## Chain 3 finished in 4.4 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 4.5 seconds.
## Total execution time: 13.8 seconds.
##
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 5.7 seconds.
## Chain 2 finished in 5.2 seconds.
## Chain 3 finished in 5.8 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 5.6 seconds.
## Total execution time: 16.9 seconds.
##
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 15.2 seconds.
## Chain 2 finished in 13.9 seconds.
## Chain 3 finished in 14.9 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 14.6 seconds.
## Total execution time: 44.2 seconds.
##
## Running MCMC with 3 sequential chains...
##
## Chain 1 finished in 7.5 seconds.
## Chain 2 finished in 7.4 seconds.
## Chain 3 finished in 7.7 seconds.
##
## All 3 chains finished successfully.
## Mean chain execution time: 7.5 seconds.
## Total execution time: 22.9 seconds.
results <- nowcasts_plot(out, models_to_run = c("fixed_q", "fixed_b", "b_poly", "b_spline"))
results$nowcasts
## [[1]]
## date case_true case_reported mean_fixed_q lower_fixed_q upper_fixed_q
## 1 2009-07-05 30 17 133.31839 46.64999 401.1522
## 2 2009-07-12 64 11 73.65665 22.21361 230.0231
## 3 2009-07-19 223 2 56.81406 11.07223 192.6796
## mean_fixed_b lower_fixed_b upper_fixed_b mean_b_poly lower_b_poly
## 1 19.802307 11.818367 42.16699 17.78680 11.230400
## 2 10.504008 4.755791 25.90012 14.33943 4.949041
## 3 7.518189 1.797386 20.99066 15.02397 2.345150
## upper_b_poly mean_b_spline lower_b_spline upper_b_spline
## 1 33.31173 17.182216 11.998743 24.95523
## 2 38.86748 8.730459 4.407380 15.99467
## 3 43.36155 5.577643 1.322419 14.54677
##
## [[2]]
## date case_true case_reported mean_fixed_q lower_fixed_q upper_fixed_q
## 1 2009-07-05 30 18 111.82049 45.33370 274.6593
## 2 2009-07-12 64 11 65.84152 24.38255 167.1970
## 3 2009-07-19 223 18 91.26846 32.52952 234.8767
## 4 2009-07-26 639 32 312.52320 100.94102 845.6314
## mean_fixed_b lower_fixed_b upper_fixed_b mean_b_poly lower_b_poly
## 1 19.98174 13.893788 31.25885 19.41739 12.739805
## 2 11.42998 6.725125 19.56273 10.07463 6.006593
## 3 15.43489 8.418845 28.31064 12.32245 7.201909
## 4 46.22332 23.913853 94.03930 31.72390 18.992755
## upper_b_poly mean_b_spline lower_b_spline upper_b_spline
## 1 43.40384 19.15419 14.016965 27.46154
## 2 20.93934 11.27083 6.629580 19.90693
## 3 25.21783 15.02600 8.223486 27.25052
## 4 69.53310 41.05433 24.780068 69.41412
##
## [[3]]
## date case_true case_reported mean_fixed_q lower_fixed_q upper_fixed_q
## 1 2009-07-05 30 19 106.16023 51.06794 231.3775
## 2 2009-07-12 64 25 89.31174 40.93426 193.8198
## 3 2009-07-19 223 91 336.06359 159.10863 732.6172
## 4 2009-07-26 639 329 2885.77450 1310.22925 6407.2708
## 5 2009-08-02 1331 272 17519.04340 7537.87350 39700.5975
## mean_fixed_b lower_fixed_b upper_fixed_b mean_b_poly lower_b_poly
## 1 312.5958 64.33369 1275.000 16487.280 1229.1370
## 2 272.9074 53.62035 1132.745 5555.680 461.5512
## 3 1015.7924 198.15353 4235.884 7117.365 672.1575
## 4 6533.0157 1269.40375 26963.682 16105.715 1621.8465
## 5 14554.9245 2734.95650 61123.365 13406.400 1292.1498
## upper_b_poly mean_b_spline lower_b_spline upper_b_spline
## 1 78793.81 22.62987 16.05667 36.6624
## 2 26405.45 43.88084 16.64506 122.0198
## 3 33576.62 295.06132 113.89650 709.6286
## 4 74073.96 1251.80798 625.86873 2495.7258
## 5 61660.34 721.40568 374.07145 1678.9457
##
## [[4]]
## date case_true case_reported mean_fixed_q lower_fixed_q upper_fixed_q
## 1 2009-07-05 30 22 101.1988 52.48704 206.0062
## 2 2009-07-12 64 39 115.6864 59.28950 235.2852
## 3 2009-07-19 223 154 532.3109 276.01545 1087.5330
## 4 2009-07-26 639 463 3296.0519 1691.39250 6701.4555
## 5 2009-08-02 1331 563 7883.8010 4012.39450 16188.9775
## 6 2009-08-09 1324 183 6701.1398 3329.85925 13867.4200
## mean_fixed_b lower_fixed_b upper_fixed_b mean_b_poly lower_b_poly
## 1 664.5866 139.2546 2720.393 20933.383 1743.4123
## 2 772.1524 161.4177 3151.671 9917.254 941.5613
## 3 3460.7446 726.0558 13839.747 17488.672 1864.3603
## 4 17821.7100 3702.5860 72474.202 36744.624 4038.6457
## 5 35965.4325 7432.5652 146923.800 31145.257 3272.3892
## 6 23525.8658 4791.7868 95778.742 8879.344 835.0984
## upper_b_poly mean_b_spline lower_b_spline upper_b_spline
## 1 112129.92 22.80778 17.03424 34.29437
## 2 52240.59 113.12357 34.57167 354.13070
## 3 91383.51 976.56980 385.76797 2345.92350
## 4 193812.32 3290.25938 1550.31900 6773.90100
## 5 157226.37 1877.17243 968.19653 4044.52550
## 6 45372.17 779.78248 269.35700 2563.51500
##
## [[5]]
## date case_true case_reported mean_fixed_q lower_fixed_q upper_fixed_q
## 1 2009-07-05 30 23 90.37367 49.29757 177.1377
## 2 2009-07-12 64 47 123.33489 66.48288 245.2817
## 3 2009-07-19 223 179 564.83845 311.26555 1107.5337
## 4 2009-07-26 639 519 2711.89706 1482.69625 5357.3300
## 5 2009-08-02 1331 749 5461.32794 2986.19050 10787.9850
## 6 2009-08-09 1324 549 5180.85360 2823.86550 10193.5200
## 7 2009-08-16 860 145 4031.76048 2121.54850 8135.5607
## mean_fixed_b lower_fixed_b upper_fixed_b mean_b_poly lower_b_poly
## 1 825.9353 144.0060 3325.223 21696.877 1252.5477
## 2 1152.8408 197.1500 4616.323 13952.930 1047.3467
## 3 5277.4571 907.5011 20748.145 29849.471 2828.3058
## 4 23824.1118 4061.5960 93118.340 65002.842 6927.0462
## 5 46710.7421 7938.0132 184770.125 62972.679 6733.1358
## 6 43012.6666 7214.7577 171637.475 29594.031 2971.3160
## 7 25582.5461 4163.9907 102242.300 9233.209 800.6328
## upper_b_poly mean_b_spline lower_b_spline upper_b_spline
## 1 123191.45 23.07897 17.81705 31.51302
## 2 80784.31 238.92735 61.40564 776.20422
## 3 163114.75 1689.78263 653.41512 3887.71025
## 4 352419.95 3467.46140 1824.79325 6798.57050
## 5 345831.37 2994.99293 1557.71750 6429.31850
## 6 163362.07 2762.66170 1295.12625 5732.20225
## 7 51634.50 918.25508 230.27995 3504.31475
##
## [[6]]
## date case_true case_reported mean_fixed_q lower_fixed_q upper_fixed_q
## 1 2009-07-05 30 25 81.99376 48.45588 148.7370
## 2 2009-07-12 64 53 125.33546 73.32550 224.3693
## 3 2009-07-19 223 198 559.78898 329.34952 1004.1915
## 4 2009-07-26 639 550 2327.84895 1367.77950 4189.2705
## 5 2009-08-02 1331 932 4500.98974 2638.69725 8079.0722
## 6 2009-08-09 1324 767 4384.64140 2566.93750 7900.2925
## 7 2009-08-16 860 417 3373.63354 1967.41375 6129.6857
## 8 2009-08-23 879 249 5833.38616 3325.97500 10632.5375
## mean_fixed_b lower_fixed_b upper_fixed_b mean_b_poly lower_b_poly
## 1 377.9449 104.0189 1499.828 12431.61 617.1515
## 2 590.4376 160.0075 2288.200 10063.84 694.3955
## 3 2667.9433 722.2312 10318.495 24016.91 2204.7238
## 4 10961.1402 2948.7370 42443.912 53547.40 6062.0782
## 5 21754.0638 5807.2428 84881.142 59098.46 7723.4870
## 6 21515.0723 5683.8283 84707.110 33273.22 4601.1780
## 7 16059.0498 4158.0995 63611.677 14523.50 1896.2160
## 8 21205.7208 5344.4838 84359.415 11451.78 1311.5397
## upper_b_poly mean_b_spline lower_b_spline upper_b_spline
## 1 67606.31 26.62506 18.72072 31.99615
## 2 54624.59 341.40807 93.16127 956.24750
## 3 126258.27 1605.09631 746.39595 3145.93750
## 4 266306.37 2292.88952 1303.80900 4530.89975
## 5 287455.47 3114.13039 2019.61825 5240.50500
## 6 165668.87 3835.33043 2008.81750 7637.41075
## 7 74280.96 2552.98415 1036.18700 5705.89825
## 8 58840.54 1297.60192 390.42020 4189.35175
results$plots
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
nowcasts_plot(out, models_to_run = c("fixed_q"))[["plots"]]
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
nowcasts_plot(out, models_to_run = c("fixed_b"))[["plots"]]
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
nowcasts_plot(out, models_to_run = c("b_poly"))[["plots"]]
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
nowcasts_plot(out, models_to_run = c("b_spline"))[["plots"]]
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
